Hands-on Exercise 3: Network Constrained Spatial Point Patterns Analysis

1 Overview

Network constrained Spatial Point Patterns Analysis (NetSPAA) is a collection of spatial point patterns analysis methods special developed for analysing spatial point event occurs on or alongside network. The spatial point event can be locations of traffic accident or childcare centre for example. The network, on the other hand can be a road network or river network.

In this hands-on exercise, we will gain hands-on experience on using appropriate functions of spNetwork package:

• to derive network kernel density estimation (NKDE), and

• to perform network G-function and k-function analysis

2 The Data

In this study, we will analyse the spatial distribution of childcare centre in Punggol planning area. For the purpose of this study, two geospatial data sets will be used. They are:

Data Sources

Type	Name	Details
Geospatial	Punggol_St	Line feature geospatial data Pertains to road network within Punggol planning area Format: SHP (ESRI Shapefile)
Geospatial	Punggol_CC	Point features geospatial data Pertains to the location of childcare centres within Punggol planning area Format: SHP (ESRI Shapefile)

3 Installing and launching the R packages

Packages

In this hands-on exercise, four R packages will be used, they are:

Package	Description
sf	Provides functions to manage, processing, and manipulate Simple Features, a formal geospatial data standard that specifies a storage and access model of spatial geometries such as points, lines, and polygons.
spNetwork	Provides functions to perform Spatial Point Patterns Analysis such as kernel density estimation (KDE) and K-function on network. It also can be used to build spatial matrices (‘listw’ objects like in spdep package) to conduct any kind of traditional spatial analysis with spatial weights based on reticular distances.
tidyverse	A collection of functions for performing data science task such as importing, tidying, wrangling data and visualising data.
tmap	For thematic mapping; provides functions for plotting cartographic quality static point patterns maps or interactive maps by using leaflet API

Code

To install and launch the four R packages.

install.packages("rgdal", repos = "https://packagemanager.posit.co/cran/2023-10-13")

pacman::p_load(sf, spNetwork, tmap, tidyverse)

4 Import Data and Preparation

4.1 Import

The code chunk below uses st_read() of sf package to important Punggol_St and Punggol_CC geospatial data sets into RStudio as sf data frames.

network <- st_read(dsn="data/geospatial", 
                   layer="Punggol_St")

Reading layer `Punggol_St' from data source 
  `C:\kytjy\ISSS626-GAA\Hands-on_Ex\Hands-on_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM

childcare <- st_read(dsn="data/geospatial",
                     layer="Punggol_CC")

Reading layer `Punggol_CC' from data source 
  `C:\kytjy\ISSS626-GAA\Hands-on_Ex\Hands-on_Ex03\data\geospatial' 
  using driver `ESRI Shapefile'
Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM

4.2 Checking Contents

Network

network

Simple feature collection with 2642 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 34038.56 ymin: 40941.11 xmax: 38882.85 ymax: 44801.27
Projected CRS: SVY21 / Singapore TM
First 10 features:
     LINK_ID                   ST_NAME                       geometry
1  116130894                PUNGGOL RD LINESTRING (36546.89 44574....
2  116130897 PONGGOL TWENTY-FOURTH AVE LINESTRING (36546.89 44574....
3  116130901   PONGGOL SEVENTEENTH AVE LINESTRING (36012.73 44154....
4  116130902   PONGGOL SEVENTEENTH AVE LINESTRING (36062.81 44197....
5  116130907           PUNGGOL CENTRAL LINESTRING (36131.85 42755....
6  116130908                PUNGGOL RD LINESTRING (36112.93 42752....
7  116130909           PUNGGOL CENTRAL LINESTRING (36127.4 42744.5...
8  116130910               PUNGGOL FLD LINESTRING (35994.98 42428....
9  116130911               PUNGGOL FLD LINESTRING (35984.97 42407....
10 116130912            EDGEFIELD PLNS LINESTRING (36200.87 42219....

str(network)

Classes 'sf' and 'data.frame':  2642 obs. of  3 variables:
 $ LINK_ID : num  1.16e+08 1.16e+08 1.16e+08 1.16e+08 1.16e+08 ...
 $ ST_NAME : chr  "PUNGGOL RD" "PONGGOL TWENTY-FOURTH AVE" "PONGGOL SEVENTEENTH AVE" "PONGGOL SEVENTEENTH AVE" ...
 $ geometry:sfc_LINESTRING of length 2642; first list element:  'XY' num [1:2, 1:2] 36547 36559 44575 44614
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA
  ..- attr(*, "names")= chr [1:2] "LINK_ID" "ST_NAME"

st_crs(network)

Coordinate Reference System:
  User input: SVY21 / Singapore TM 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

Childcare

childcare

Simple feature collection with 61 features and 1 field
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 34423.98 ymin: 41503.6 xmax: 37619.47 ymax: 44685.77
z_range:       zmin: 0 zmax: 0
Projected CRS: SVY21 / Singapore TM
First 10 features:
      Name                      geometry
1   kml_10 POINT Z (36173.81 42550.33 0)
2   kml_99 POINT Z (36479.56 42405.21 0)
3  kml_100 POINT Z (36618.72 41989.13 0)
4  kml_101 POINT Z (36285.37 42261.42 0)
5  kml_122  POINT Z (35414.54 42625.1 0)
6  kml_161 POINT Z (36545.16 42580.09 0)
7  kml_172 POINT Z (35289.44 44083.57 0)
8  kml_188 POINT Z (36520.56 42844.74 0)
9  kml_205  POINT Z (36924.01 41503.6 0)
10 kml_222 POINT Z (37141.76 42326.36 0)

str(childcare)

Classes 'sf' and 'data.frame':  61 obs. of  2 variables:
 $ Name    : chr  "kml_10" "kml_99" "kml_100" "kml_101" ...
 $ geometry:sfc_POINT of length 61; first list element:  'XYZ' num  36174 42550 0
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
  ..- attr(*, "names")= chr "Name"

st_crs(childcare)

Coordinate Reference System:
  User input: SVY21 / Singapore TM 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

5 Visualising Geospatial data

Using Plot

par(bg = '#E4D5C9')

plot(st_geometry(network))

plot(childcare,
     add=T,
     col='#800200',
     pch = 19)

par(bg = '#E4D5C9')

plot(network)

plot(childcare,
     add=T,
     col='#800200',
     pch = 19)

Using Tmap

For achieving a visually appealing and interactive representation of geospatial data, the tmap package’s mapping function can be employed.

tmap_mode('view')

tm_shape(childcare) + 
  tm_dots(col = "#800200") + 
  tm_shape(network) +
  tm_lines()

6 Network Constrained KDE (NetKDE) Analysis provided in spNetwork

6.1 Preparing the lixels objects

Prior to computing NetKDE, it is necessary to partition the SpatialLines object into lixels with a specified minimum distance. This operation can be accomplished using the lixelize_lines() function from the spNetwork package.

Lixel Length 700

lixels <- lixelize_lines(lines = network, #<< SpatialLinesDataFrame
                         lx_length = 700, #<< Length of a lixel
                         mindist = 350) #<< Minimum length of a lixel

Lixel Length 750

lixels_750 <- lixelize_lines(lines = network, #<< SpatialLinesDataFrame
                         lx_length = 750, #<< Length of a lixel
                         mindist = 375) #<< Minimum length of a lixel

Lixelize_lines Function

• Dimensions for Lixels Objects:

  • Set the length of a lixel (lx_length) to 700m.

  • Set the minimum length of a lixel (mindist) to 350m.

  • After cut, if the final lixel is shorter than the minimum distance, it will be added to the previous lixel.

  • Segments that are already shorter than the minimum length are not modified.

  • If the minimum length is NULL, then mindist = maxdist/10.

  • Additional Information about Lixelize_lines Function:

  • Lixelize_lines is used to cut a SpatialLines object into lixels with a specified minimal distance.

  • The function also supports multicore processing through lixelize_lines.mc().

• Post-cut Considerations:
  • After cutting, if the length of the final lixel is shorter than the minimum distance, it is added to the previous lixel.
  • If the minimum distance is NULL, then mindist is set to maxdist/10.
  • Segments that are already shorter than the minimum distance are not modified.

6.2 Generate line centre points using lines_center() of spNetwork

• lines_center() of spNetwork is used to generate a SpatialPointsDataFrame with line center points.

• Points are located at center of line based on the length of the line.

Lixel Length 700

samples <- lines_center(lixels)

Lixel Length 750

samples_750 <- lines_center(lixels_750)

6.3 Performing NetKDE

To compute the NetKDE:

Lixel Length 700

densities <- nkde(network, 
                  events = childcare,
                  w = rep(1, nrow(childcare)),
                  samples = samples,
                  kernel_name = "quartic", #<<
                  bw = 300, 
                  div= "bw", 
                  method = "simple", #<< Can be simple, discontinuous, continuous
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  #agg = 5, #<< Aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

Lixel Length 750

densities_750 <- nkde(network, 
                  events = childcare,
                  w = rep(1, nrow(childcare)),
                  samples = samples_750,
                  kernel_name = "quartic", #<<
                  bw = 300, 
                  div= "bw", 
                  method = "simple", #<< Can be simple, discontinuous, continuous
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  #agg = 5, #<< Aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

Kernel Method and Arguments:

• The code chunk reveals the use of the quartic kernel (kernel_name argument).

• spNetwork supports various kernel methods, including triangle, gaussian, scaled gaussian, tricube, cosine, triweight, epanechnikov, or uniform.

Calculation Methods for NKDE:

• The method argument indicates the use of the “simple” method for calculating NetKDE.

• spNetwork offers three methods for NKDE:

  • simple: Distances between events and sampling points are replaced by network distances. The kernel formula is adjusted to calculate density over a linear unit instead of an areal unit.

  • discontinuous: Proposed by Okabe et al. (2008), this method “divides” the mass density of an event at intersections of lixels.

  • continuous: An alternative version proposed by Okabe et al. (2008) adjusts the density before intersections, making the function continuous.

User Guide Reference:

• The user guide of the spNetwork package provides a comprehensive discussion of nkde(). It is recommended to read the guide to understand various parameters for calibrating the NetKDE model.

Additional Notes on Arguments:

• The chosen kernel method is quartic, and the decision is explained.

• spNetwork supports alternative kernel methods such as triangle, gaussian, scaled gaussian, tricube, cosine, triweight, epanechnikov, or uniform.

• The selected method for NKDE calculation is “simple,” and the reasons for its use are explained.

• Other supported methods include “discontinuous” and “continuous,” each with specific characteristics described in the code chunk.

6.4 Visualising NetKDE

6.4.1 Insert computed density values (i.e. densities) into samples and lixels objects as density field

Lixel Length 700

samples$density <- densities
lixels$density <- densities

Lixel Length 750

samples$density_750 <- densities_750
lixels$density_750 <- densities_750

6.4.2 Rescale density values from number of events per m to number of events per km

As the svy21 projection system is in meters, the resulting density values are very small (e.g., 0.0000005). The code below employed to rescale the density values from the number of events per meter to the number of events per km.

Lixel Length 700

# rescaling to help the mapping

samples$density <- samples$density*1000
lixels$density <- lixels$density*1000

Lixel Length 750

samples$density_750 <- samples$density_750*1000
lixels$density_750 <- lixels$density_750*1000

`

6.4.3 Using tmap package to plot map after rescaling

tmap packages can be used to prepare interactive and high cartographic quality map visualisation.

Lixel Length 700

tmap_mode('view')
tm_shape(lixels)+
  tm_lines(col="density")+
tm_shape(childcare)+
  tm_dots()

Lixel Length 750

tm_shape(lixels)+
  tm_lines(col="density_750")+
tm_shape(childcare)+
  tm_dots()

Interpretation

• Road segments with relatively higher density of childcare centres (darker color)
• Road segments with relatively lower density of childcare centres (lighter color)

7 Network Constrained G- and K-Function Analysis

Complete Spatial Randomness (CSR) Test

• Objective: Conducting CSR test using the kfunctions() function from the spNetwork package.

• Null Hypothesis (\(H_0\)): The observed spatial point events (i.e., distribution of childcare centres) exhibit a uniform distribution over a street network in Punggol Planning Area.

• CSR Test Assumption:

  • The CSR test relies on the assumption of a binomial point process.

  • Assumption implies that childcare centres are randomly and independently distributed over the street network.

• Interpretation of Results:

  • If the null hypothesis is rejected:

    • Inference: The distribution of childcare centres shows spatial interactions and dependence.

    • Resultant Patterns: Nonrandom patterns may be observed.

• CSR Test Execution:

  • Execution involves utilizing the kfunctions() function from the spNetwork package.

Code

kfun_childcare <- kfunctions(network, 
                             childcare,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 50, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05)

Insights

Arguments Used:

• Ten arguments are employed in the code chunk, namely:

  • lines: A SpatialLinesDataFrame with sampling points.
  • points: A SpatialPointsDataFrame representing points on the network.
  • start: Start value for evaluating the k and g functions.
  • end: Last value for evaluating the k and g functions.
  • step: Jump between two evaluations of the k and g functions.
  • width: Width of each donut for the g-function.
  • nsim: Number of Monte Carlo simulations (50 simulations in the example).
  • resolution: Resolution for simulating random points on the network.
  • conf_int: Width of the confidence interval (default = 0.05).
  • For additional arguments, refer to the user guide of the spNetwork package.

Output of kfunctions():

The function outputs a list with the following components:

-   plotkA: ggplot2 object representing k-function values.
-   plotgA: ggplot2 object representing g-function values.
-   valuesA: DataFrame with values used to build the plots.

We can visualise the ggplot2 object of k-function by using the code chunk below.

kfun_childcare$plotk + 
  labs(title ="K-Function") +
  theme(panel.background = element_rect(fill = "#E4D5C9"),
        plot.background = element_rect(fill = "#E4D5C9"),
        panel.grid.major = element_line(colour = "#efe7df", linetype = 1, linewidth = 0.5),
        panel.grid.minor = element_line(colour = "#efe7df", linetype = 1, linewidth= 0.5),
        plot.title = element_text(face = "bold", size = 12, hjust = 0.5),
        )

Observations from the Graph:

• The blue line indicates the empirical network K-function for childcare centers in Punggol.
• A gray envelope represents results from 50 simulations spanning the 2.5% to 97.5% interval.

Inference:

• Blue line values between 125m-687.5m fall below the gray envelope.

Conclusion:

• Implies that childcare centers in Punggol exhibit a regular pattern within the 125m-687.5m distance range.

Kfunctions values:

8 Reference

Kam, T. S. Network Constrained Spatial Point Patterns. R for Geospatial Data Science and Analytics. https://r4gdsa.netlify.app/chap07.html